This notebook is part of PyBroMo a python-based single-molecule Brownian motion diffusion simulator that simulates confocal smFRET experiments.
In this notebook we generate single-polulation static smFRET data files from the same diffusion trajectories. These files are needed by the next notebook to generate smFRET data with 2-state dynamics.
In [ ]:
%matplotlib inline
from pathlib import Path
import numpy as np
import tables
import matplotlib.pyplot as plt
import seaborn as sns
import pybromo as pbm
import phconvert as phc
print('Numpy version:', np.__version__)
print('PyTables version:', tables.__version__)
print('PyBroMo version:', pbm.__version__)
print('phconvert version:', phc.__version__)
In [ ]:
SIM_PATH = 'data/'
We assume a $\gamma = 0.7$ and two populations, one with $E_{PR}=0.75$ and the other $E_{PR}=0.4$
In [ ]:
Epr1 = 0.75
Epr2 = 0.4
The corrected $E$ for the two populations are:
In [ ]:
gamma = 0.7
E1 = Epr1 /(Epr1 * (1 - gamma) + gamma)
E2 = Epr2 /(Epr2 * (1 - gamma) + gamma)
E1, E2
The simulation takes the uncorrected $E_{PR}$ as input.
We want to simulate a second population that has measured brightness scaling as if the difference was only due to the $\gamma$ factor. Using the definitions $\Lambda$ and $\Lambda_\gamma$ from (Ingargiola 2017), we can use the relation:
$$\frac{E}{E_{PR}} = \frac{\Lambda}{\Lambda_\gamma}$$Solving for $\Lambda_\gamma$ or $\Lambda$ we get:
$$ \Lambda_\gamma = \Lambda\frac{E_{PR}}{E}$$$${\Lambda} = {\Lambda_\gamma}\frac{E}{E_{PR}} $$Since $\Lambda_\gamma$ is gamma corrected does not depend on $E$. We can compute it from the parameters of population 1 and then use it for finding $\Lambda$ for population 2:
In [ ]:
Λ1 = 200e3 # kcps, the detected (i.e. uncorrected) peak emission rate for population 1
Λγ = Λ1 * Epr1 / E1
Λ2 = Λγ * E2 / Epr2
Λ2
In [ ]:
Λ2 = np.round(Λ2, -3)
Λ2
In [ ]:
S = pbm.ParticlesSimulation.from_datafile('0eb9', mode='a', path=SIM_PATH)
In [ ]:
S.particles.diffusion_coeff_counts
In [ ]:
#S = pbm.ParticlesSimulation.from_datafile('44dc', mode='a', path=SIM_PATH)
In [ ]:
params1 = dict(
em_rates = (Λ1,), # Peak emission rates (cps) for each population (D+A)
E_values = (Epr1,), # FRET efficiency for each population
num_particles = (35,), # Number of particles in each population
bg_rate_d = 900, # Poisson background rate (cps) Donor channel
bg_rate_a = 600, # Poisson background rate (cps) Acceptor channel
)
We can now define population 2:
In [ ]:
params2 = dict(
em_rates = (Λ2,), # Peak emission rates (cps) for each population (D+A)
E_values = (Epr2,), # FRET efficiency for each population
num_particles = (35,), # Number of particles in each population
bg_rate_d = 900, # Poisson background rate (cps) Donor channel
bg_rate_a = 600, # Poisson background rate (cps) Acceptor channel
)
Finally, we also define a static mixture of the two populations:
In [ ]:
params_mix = dict(
em_rates = (Λ1, Λ2), # Peak emission rates (cps) for each population (D+A)
E_values = (Epr1, Epr2), # FRET efficiency for each population
num_particles = (20, 15), # Number of particles in each population
bg_rate_d = 900, # Poisson background rate (cps) Donor channel
bg_rate_a = 600, # Poisson background rate (cps) Acceptor channel
)
Population 1: Create the object that will run the simulation and print a summary:
In [ ]:
mix_sim = pbm.TimestapSimulation(S, **params1)
mix_sim.summarize()
Run the simualtion:
In [ ]:
rs = np.random.RandomState(1234)
mix_sim.run(rs=rs, overwrite=False, skip_existing=True)
Save simulation to a smFRET Photon-HDF5 file:
In [ ]:
mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola',
author_affiliation='UCLA'))
Population 2: Create the object that will run the simulation and print a summary:
In [ ]:
mix_sim = pbm.TimestapSimulation(S, **params2)
mix_sim.summarize()
In [ ]:
rs = np.random.RandomState(1234)
mix_sim.run(rs=rs, overwrite=False, skip_existing=True)
In [ ]:
mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola',
author_affiliation='UCLA'))
Static mixture: Create the object that will run the simulation and print a summary:
In [ ]:
mix_sim = pbm.TimestapSimulation(S, **params_mix)
mix_sim.summarize()
In [ ]:
rs = np.random.RandomState(1234)
mix_sim.run(rs=rs, overwrite=False, skip_existing=True)
In [ ]:
mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola',
author_affiliation='UCLA'))
In [ ]:
!rsync -av --exclude 'pybromo_*.hdf5' /mnt/archive/Antonio/pybromo /mnt/wAntonio/dd
The generated Photon-HDF5 files can be analyzed by any smFRET burst analysis program. Here we show an example using the opensource FRETBursts program:
In [ ]:
import fretbursts as fb
In [ ]:
filepath = list(Path(SIM_PATH).glob('smFRET_*'))
filepath
In [ ]:
d = fb.loader.photon_hdf5(str(filepath[0]))
In [ ]:
d
In [ ]:
d.A_em
In [ ]:
fb.dplot(d, fb.timetrace);
In [ ]:
d.calc_bg(fun=fb.bg.exp_fit, tail_min_us='auto', F_bg=1.7, time_s=5)
In [ ]:
fb.dplot(d, fb.timetrace_bg)
In [ ]:
d.burst_search(F=7)
In [ ]:
d.num_bursts
In [ ]:
ds = d.select_bursts(fb.select_bursts.size, th1=20)
In [ ]:
ds.num_bursts
In [ ]:
with plt.rc_context({#'font.size': 10,
#'savefig.dpi': 200,
'figure.dpi': 150}):
for i in range(3):
fig, ax = plt.subplots(figsize=(100, 3))
fb.dplot(d, fb.timetrace, binwidth=0.5e-3, tmin=i*10, tmax=(i+1)*10, bursts=True,
plot_style=dict(lw=1), ax=ax);
ax.set_xlim(i*10, (i+1)*10);
display(fig)
plt.close(fig)
In [ ]:
fb.dplot(ds, fb.hist_fret, pdf=False)
plt.axvline(0.4, color='k', ls='--');
In [ ]:
fb.bext.burst_data(ds)
In [ ]:
fb.dplot(d, fb.hist_size)
In [ ]:
fb.dplot(d, fb.hist_width)